clc ;clear
clearvars
clear global
close all

tic

%% set parameters

% structural parameters
alp   = 1.0;            % input suitability (Changed from 1)
m_phi = 0.;          % mean of fundamental productivity
m_del = 0.;          % mean of fundamental quality
v_phi = 1.;          % std of fundamental productivity
v_del = 1.;          % std of fundamental quality
v_cor = 0;           % correlation between phi and del

s_psi = 0;           % elasticity of psi function
s_xi  = .9565;       % measure of xi variance (shape parameter of the weibull distribution)
L     = 1;           % endowment of emission-intensive inputs

m_psi = 0.216;
m_psi_II = m_psi;      % scale of psi function (average relationship cost)
m_psi_IC = m_psi;      % average relationship cost of firms in different region (domestic)
% m_psi_IFI = 0.216;     % average relationship cost of inland firms between home and foreign countries
% m_psi_IFC = 0.216;     % average relationship cost of inland and coastal firms between home and foreign countries
% m_psi_CFC = 0.216;     % average relationship cost of coastal firms between home and foreign countries
m_psi_F = 0.6;

sig = 4;            % elasticity of substitution

%% Key Parameters
a_LC   = 0.06; %0.11; %         % average relocation cost (changed from 0.02 to 0.1 to 0.25)
tau_B  = 1.2;           % international trade cost
tau_D  = 1.17; %1.17;         %1.09     % domestic trade cost in home country (1.1 work)
tau_DF = 1.17; %1.17;        % domestic trade cost in foreign country
w      = 1.0;           % unit price of emission-intensive inputs
t      = 0.134; % 0.439;       % emission tax rate 
eps    = 0.000094; % 0.034;        % scale factor of emission-intensive inputs on emissions (Carbon content in natural gas for example)
net_open = 1;        % network openness
%%%%%%%%%%%%%%%%%%

M = [m_phi;m_del];                                                  % mean vector of (log(phi),log(del))
V = [v_phi^2,v_cor*v_phi*v_del;                                     % covariance matrix of (log(phi),log(del))
     v_cor*v_phi*v_del,v_del^2]; 
g_chi_fun = @(phi,del) mvnpdf([log(phi);log(del)],M,V)/(phi*del);   % density function of (phi,del)

% solver parameters
qt_min = .1;        % min quantile for phi/del distribution
qt_max = .9;        % max quantile for phi/del distribution
n_grid = 15;   % grid size
tol    = 1E-3;      % numerical tolerance for iterations
max_iter = 200;     % max iterations


%% store parameters

% structural parmaeters
struc_params.sig   = sig;
struc_params.alp   = alp;
struc_params.m_phi = m_phi;
struc_params.m_del = m_del;
struc_params.v_phi = v_phi;
struc_params.v_del = v_del;
struc_params.v_cor = v_cor;
struc_params.m_psi = m_psi;
struc_params.m_psi_II = m_psi_II;
struc_params.m_psi_IC = m_psi_IC;
struc_params.m_psi_F = m_psi_F;
struc_params.s_psi = s_psi;
struc_params.s_xi  = s_xi;
struc_params.L     = L;
struc_params.tau_B = tau_B;
struc_params.tau_D = tau_D;
struc_params.tau_DF = tau_DF;
struc_params.w     = w;
struc_params.t     = t;
struc_params.eps   = eps;
struc_params.a_LC  = a_LC;
struc_params.net_open  = net_open;


% solver parameters
solv_params.qt_min   = qt_min;
solv_params.qt_max   = qt_max;
solv_params.n_grid   = n_grid;
solv_params.tol      = tol;
solv_params.max_iter = max_iter;

TAUS = [1.0, 1.025, 1.05, 1.075, 1.1, 1.125, 1.15, 1.175, 1.2];

OUTPUT = [];
solv_params.matching_flag = 'endogenous';
high_trade_cost = solver_exploc_matrix(struc_params, solv_params, 1.7);

avg_MC = []; 
avg_MCF = [];
mean_matching_prob = [];
entry_cost = zeros(n_grid, n_grid, length(TAUS));

endog_results = [];


parfor (idx = 1:length(TAUS), length(TAUS))
    disp("parallel loop initialized")
    tau = TAUS(idx);
    %% solve model and load outcomes
    disp("solve model")
    model_exploc1_matrix9 = solver_exploc_matrix(struc_params, solv_params, tau, 0);
    disp("save results")
    
    endog_results = [endog_results, model_exploc1_matrix9];

    tmp = model_exploc1_matrix9.RC_total;
    tmp(tmp == 0) = min(min(tmp(tmp > 0)));
    entry_cost(:, :, idx) = tmp;
    avg_MCF = [avg_MCF, model_exploc1_matrix9.avg_m_C_F];
    avg_MC = [avg_MC, model_exploc1_matrix9.avg_m_C_d];
    mean_matching_prob = [mean_matching_prob, model_exploc1_matrix9.mean_matchingprob];

    out = excel_counterfactual(model_exploc1_matrix9, tau)
    disp('concat results')
    OUTPUT = [OUTPUT,out];
end

disp("reshape")
OUTPUT = OUTPUT';

%% Run the IO Counterfactual
disp("get mean matching probability")
avg_MCF = mean(avg_MCF);
avg_MC = mean(avg_MC);
mean_matching_prob = mean(mean_matching_prob);

struc_params.io_param = mean_matching_prob;
struc_params.io_param_D = avg_MC;
struc_params.io_param_F = avg_MCF / avg_MC;
struc_params.entry_cost = entry_cost;
solv_params.matching_flag = 'exogenous';

IO_CF = [];
exog_results = [];


parfor (idx = 1:length(TAUS), length(TAUS))
    tau = TAUS(idx);
    entry_cost_exog = entry_cost(:, :, idx);
    %% solve model and load outcomes
    % solve model
    entry_status_exog = endog_results(idx).entry_status; % Fix the set of firms in the experiment. 
    model_exploc1_matrix9 = solver_exploc_matrix(struc_params, solv_params, tau, entry_status_exog); 
    
    exog_results= [exog_results, model_exploc1_matrix9];

    IO_CF = [IO_CF, excel_counterfactual(model_exploc1_matrix9, tau)];
end
IO_CF = IO_CF';

NAMES = ["Simulation Key", "tau", "mean_matchingprob", "share_coastal", "share_inland", "inland_eint",...
            "coastal_eint", "share_inland_etotal", "share_coastal_etotal", "total_emission",...
            "total_eint", "share_inland_domestic", "share_inland_exporter",...
            "share_coastal_domestic", "share_coastal_exporter", "share_dom",...
            "share_intl", "eint_ratio_coastal_inland", "eint_ratio_inland_coastal",...
            "N_coastal", "N_inland", "avg_m_S_d", "avg_m_S_F", "avg_m_C_F", "avg_m_C_d",...
            "inland_etotal", "coastal_etotal"];


OUTPUT = [repelem("Endog", length(TAUS))', OUTPUT];
IO_CF = [repelem("IO", length(TAUS))', IO_CF];


% Save to Excel
filename = 'counterfactual.xlsx';
writematrix([NAMES; OUTPUT; IO_CF],filename,'Sheet',1)%; NO_NW_CF

save 'counterfactual_results.mat'

toc